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Abstract 


>This  memorandum  is  the  second  of  three  which  describe  a  series  of  computer 
programs  for  the  calculation  of  the  growth  of  a  boundary  layer  on  a  ship  hull.  The  first 
described  computer  programs  for  the  numerical  approximation  of  ship  hulls  using 
smoothing  splines.  This  memorandum  describes  the  calculation  of  the  potential  flow 
around  the  numerically  approximated  hulls.  The  third  memorandum  will  describe  the 
calculation  of  the  boundary  layer  growth. 

The  potential  flow  is  calculated  using  a  first-order  panel  method  which  yields  the 
cartesian  components  of  the  flow  vector  at  a  number  of  control  points  on  the  hull.  These 
discrete  cartesian  flow  vectors  are  converted  to  contravariant  velocity  components  and 
splined  to  yield  a  continuous  representation  of  the  potential  flow  over  the  hull.  Zero 
Froude  number  is  assumed. 

User  guides  are  given  for  programs  implementing  the  panel  method  and  for  a  program 
which  splines  the  calculated  velocities.  A  c^c.:  (  ' 


Resume 


Le  present  memo ire  est  le  deuxieme  de  trois  memo ires  qui 
decrivent  une  sdrie  de  programmes  informat iques  permettant  de  calcuier  la 
croissance  d'une  couche  limite  sur  une  coque  de  bateau.  Le  premier 
decrivait  des  programmes  inf ormatiques  d*  approximation  numerique  des 
coques  de  bateau  &  partir  de  fonctions  splines  de  lissage.  Ce  deuxieme 
porte  sur  le  calcul  de  l’ecoulement  potential  autour  des  coques  aprochdes 
numeriquement .  Le  troisieme  portera  sur  le  calcul  de  la  croissance  des 
couches  limites. 

L'ecoulement  potential  est  etabli  au  moyen  d'une  mithode  de 
panneaux  de  premier  ordre  qui  donne  les  composantes  cartisiennes  du 
vecteur  d'ecoulement  an  certains  points  de  contrdle  de  la  coque.  Ces 
vecteurs  d'ecoulement  cartesians  discrets  sont  transforms  an  composantes 
de  vitesse,  lesquelles  sont  lissdas  par  fonctions  splines  an  une 
representation  continue  de  l'ecoulement  potentiel  autour  de  la  coque.  On 
suppose  un  nombre  de  Froude  nul. 

Des  guides  de  1* utilisation  sont  distribues;  ils  portent  sur  les 
progrananes  d* application  de  la  methods  des  panneaux  et  sur  un  programme 
de  lissage  par  splines  des  vitesses  calculdes. 
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NOMENCLATURE 


element  of  influence  coefficient  matrix 
the  p**  B-spline  of  order  k 

hull  surface  vector,  basis  for  contravariant  vectors 

hull  surface  vector,  basis  for  covariant  vectors 
orders  of  splines  in  the  x  and  s  directions  respectively 
number  of  B-splines  in  the  x  and  s  directions  respectively 
freestream  Mach  number 

normal  to  a  boundary  or  panel 
panel  control  point 

corner  points  of  non-planar  panel  (i  *  1  ,...4) 
corner  points  of  planar  panel  (i  *  1  ,...4) 

Planar  panel  definition  point 
displacement  vector  between  two  points 
distance  between  two  points.  The  magnitude  of  r. 
boundary  of  the  domain  ft 

knot  sequences  for  splines  in  x  and  s  directions  respectively 

solid  angle  normalized  by  47T 

vectors  joining  the  corners  of  the  planar  panel 

local  velocity  vector 

free  stream  velocity 

contravariant  component  of  V 
non-orthogonal  surface  coordinate  system 


x',y\z’  -  cartesian  coordinate  system  used  in  EN967 
X,Y,2  -  cartesian  coordinate  system  used  in  SMHULL 

“pq^pq  "  B-spline  coefficients  for  Y  and  Z  respectively 

i\kpq  -  B-spline  coefficients  for  Vk 

pkq  -  B-spline  coefficient  defined  in  equation  (5.7) 

cr  -  source  density 

$  -  velocity  potential 

ft  -  flow  region,  or  domain 

ABBREVIATIONS 

AP  -  aft  perpendicular 

FP  -  forward  perpendicular 

LBP  -  length  between  perpendiculars 

BOLD  CHARACTERS 

Bold  characters  are  reserved  for  vectors. 

SUBSCRIPTS 

Subscripts  i  and  j  are  generally  associated  with  the  discretized  flow  model.  Thus,  V, 
is  the  local  velocity  at  the  ith  control  point  NOT  the  itt'  component  of  the  velocity  vector. 

Subscripts  k  and  n  are  generally  associated  with  the  covariant  and  contravariant 
coordinates  of  vectors  and  the  basis  vectors  for  the  hull  coordinate  system.  Thus  vk  is 
the  kth  contravariant  component  of  vector  V  and  ek  is  the  k*1  contravariant  basis  vector. 


1  INTRODUCTION 


This  memorandum  is  the  second  of  three  which  present  programs  leading  to  the 
prediction  of  boundary  layers  starting  from  hull  section  offsets.  The  purpose  of  the 
programs  EN967  and  VELSPL  described  here  is  to  calculate  the  velocity  field  of  a  steady, 
incompressible,  potential  flow  around  a  ship  hull  whose  geometry  has  been  specified  by 

the  program  SMHULL1.  This  flow  field  may  then  be  used  as  input  to  the  program  BLAYER2 
which  calculates  the  boundary  layer  growth  on  the  hull.  Figure  illustrates  how  the 
programs  described  here  fit  into  the  overall  boundary  layer  calculation.  Since  input  data 
for  the  programs  described  here  must  be  pre-processed  using  the  programs  described  in 

the  first  report  published1,  a  reasonable  level  of  familiarity  with  that  document  is 
assumed. 

Fluid  flow  is  described  as  a  potential  flow  when  the  velocity  field  is  the  gradient  of 
a  scalar  potential.  Such  a  flow  is  irrotational  since  the  vorticity  vanishes  everywhere.  If 
the  flow  is  also  incompressible,  the  scalar  potential  is  governed  by  Laplace's  equation 
and  may  be  determined  uniquely  if  appropriate  boundary  conditions  are  specified.  The 
boundary  conditions  appropriate  for  the  flow  past  a  ship  moving  at  constant  speed  are 
that  the  flow  at  infinity  be  uniform  and  that  there  is  no  flux  of  fluid  through  the  hull  of  the 
ship.  An  additional  simplification  is  made  by  assuming  the  Froude  number  to  be  zero;  that 
is,  there  is  no  wavemaking,  and  the  free  surface  is  replaced  by  a  rigid,  reflecting  surface. 
This  approximation  is  satisfactory  except  very  near  the  surface  and  considerably 
simplifies  the  subsequent  boundary  layer  calculation.  The  flow  thus  calculated  gives  a 
good  approximation  to  the  flow  around  a  ship  hull  outside  the  boundary  layer,  and,  where 
the  boundary  layer  is  considered  'thin'  (a  reasonable  assumption  over  the  greater  part  of 
a  ship  hull),  may  be  used  to  predict  dynamic  pressures  and  streamlines  with  acceptable 
accuracy.  It  is  used  in  the  present  instance  to  provide  surface  velocities  for  a  boundary 
layer  prediction  method. 

The  calculation  of  the  incompressible  potential  flow  around  a  hull  can  be  broken  into 
two  main  steps.  The  flow  is  first  calculated  using  a  first  order  panel  method  which  yields 
values  of  the  velocity  vector  at  prescribed  control  points  on  the  hull.  The  velocity 
vectors  are  given  in  cartesian  coordinates.  When  used  for  subsequent  calculations  this 
discrete  cartesian  representation  of  the  potential  flow  is  often  not  convenient,  for  three 
reasons. 

1)  It  may  be  necessary  to  know  the  potential  flow  at  points  other  than  the 
control  points.  Some  means  of  interpolation  is  then  required. 

2)  Since  the  hull  form  is  known,  only  two  components  of  the  flow  vector  are 
needed  to  determine  it  uniquely.  Hence  the  flow  vectors  are 
overdetermined,  since  the  cartesian  representation  uses  three 
components  to  represent  the  flow  vectors. 

3)  The  cartesian  coordinate  system  is  not  often  used  in  calculations 
involving  the  potential  flow.  Transformations  between  different  coordinate 
systems  are  then  necessary. 

To  avoid  these  difficulties,  the  second  calculation  step  is  invoked;  the  flow  vectors 
at  the  control  points  are  represented  in  terms  of  a  hull-surface  coordinate  system  which 
requires  only  two  components  to  determine  the  flow  vector,  and  a  smooth  interpolation  of 
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the  flow  vectors  is  found  using  B-spline  basis  functions.  The  resulting  representation  of 
the  flow  is  much  more  convenient  to  use  in  subsequent  calculations,  and  yields  a  value 
for  the  potential  flow  vector  at  any  point  on  the  hull. 

The  calculation  of  the  discrete  cartesian  representation  of  the  potential  flow  may 
itself  be  broken  into  several  steps  each  of  which  is  done  by  a  single  program  in  the  suite 
of  programs  called  EN967. 

1)  The  hull  is  first  approximated  by  a  large  number  (usually  about  1000)  of 
planar  quadrilateral  panels.  The  program  EN967S  takes  hull  data  in  the 

form  of  B-spline  coefficients  generated  by  program  SMHULL1  and 
calculates  corner  points  and  the  control  points  (near  the  centroids)  of  the 
panels. 

2)  Each  panel  is  considered  to  be  covered  by  a  uniform  Rankine  source 
distribution  of  unit  strength.  Each  panel  thus  causes  an  incompressible 
potential  flow  which  vanishes  at  infinity.  Each  panel  also  causes  a  flux  of 
fluid  through  each  other  panel.  The  program  EN967C  calculates  a  matrix  of 
"influence  coefficients"  a^  which  give  the  flux  through  panel  i  due  to  the 
source  distribution  on  panel  j. 

3)  Each  panel  is  assigned  an  undetermined  source  strength,  O',.  The  solution 
to  the  potential  flow  problem  can'  then  be  reduced  to  the  solution  of  a 
linear  system  of  equations  to  determine  values  for  the  source  strengths 
such  that  the  no-flux  boundary  condition  is  satisfied.  The  solution  of  the 
linear  system  is  performed  in  EN967D. 

4)  When  the  panel  source  strengths  are  known,  the  velocity  at  any  point  in 
the  fluid  may  be  determined.  In  EN967E  the  velocity  (in  cartesian 
coordinates)  is  determined  at  the  control  point  of  each  panel. 

A  user's  guide  for  EN967  is  contained  in  Appendix  A.  The  EN967  suite  of  programs 
contains  additional  programs  (EN967A,  EN967B,  etc.)  which  are  used  when  calculating 
potential  flows  around  more  complex  geometries;  only  those  programs  which  are  needed 
for  the  ship  hull  problem  are  discussed  in  this  memorandum. 

The  velocity  vectors  calculated  in  EN967E  are  splined  in  hull  surface  coordinates  by 
program  VELSPL,  a  users's  guide  for  which  is  contained  in  Appendix  B. 

It  will  be  noted  in  Figure  1  that  the  present  implementation  of  the  boundary  layer 
calculation  consists  of  a  single  pass  through  the  programs.  This  is,  in  some  respects,  a 
'first  order'  calculation  since  the  effect  of  the  boundary  layer  on  the  potential  flow 
velocities  is  not  accounted  for.  Further  developments  will  incorporate  boundary  layer 
corrections  in  an  iterative  scheme  and  will  also  require  extension  of  the  present  thin 
(with  respect  to  the  principal  radii  of  curvature  of  the  surface)  boundary  layer  model  to 
thick  attached  boundary  layers  and  possibly  separated  flow  in  the  region  of  the  stern. 
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2  HULL  REPRESENTATION  AND  COORDINATE  SYSTEM 


Different  cartesian  coordinate  systems  are  used  for  the  hull  representation  in 

program  SMHULL1  and  in  the  EN967  programs,  in  addition,  a  non-orthogonal  hull  surface 
coordinate  system  is  used  by  SMHULL  and  by  program  VELSPL  when  splining  the 
calculated  velocities. 


SMHULL  generates  a  functional  relationship  between  the  cartesian  hull  coordinate 
system  (X,Y,Z)  and  the  non-orthogonal  surface  coordinate  system  (x,s).  These  two 
coordinate  systems  are  shown  in  Figure  2.  The  cartesian  system  consists  of  X  directed 
aft,  Y  to  starboard  and  Z  upwards,  with  the  origin  on  the  keel-line  at  the  FP.  The  surface 
coordinate  system  consists  of  x,  station  number,  from  0  at  the  FP  to  20  at  the  AP,  and  s, 
non-dimensionalized  arclength  around  the  girth,  from  0.0  at  the  deck  or  waterline  to  1.0 
at  the  keel. 


X  is  simply  given  by 


X(x)  = 


x  LBP 
20 


(2.1) 


and  the  other  cartesian  coordinates  are  related  to  the  surface  coordinates  by 


Nx  Nt 

Y(x,s)  =  EE  aWBP,kx<X>Bq.tc,(s)  > 


p=1  q»1 


(2.2) 


Nx  N, 

Z(x,s)  =  EE  ^Ws) 


P=1  q*1 


(2.3) 


where  Bp  k  is  the  pth  B-spline  of  order  k  and  aM  and  (3^  are  the  B-spline  coefficients  for  Y 
an'1  Z  respectively.  Details  of  this  procedure  are  given  by  Hally1  and  the  theory  and 
application  of  B-splines  are  discussed  by  de  Boor3  and  by  Hally4. 

It  is  important  to  note  that  EN967S  assumes  that  the  hull  has  been  defined  only  so 
far  as  the  waterplane  and  will  take  the  waterplane  to  be  the  maximum  value  of  Z  at  each 
section.  Thus,  hull  data  which  has  been  taken  further,  to  the  deck  for  example,  must  be 
truncated  at  the  waterplane  before  processing  by  SMHULL. 

EN967  uses  the  cartesian  coordinate  system  (x'.y'.z')  shown  in  Figure  3.  x'  is 
directed  aft,  y'  downwards  and  z'  to  starboard  with  the  origin  at  the  FP  on  the 
waterplane.  Furthermore,  as  shown  in  the  figure,  two  planes  of  symmetry  are  present: 
the  centerplane,  because  of  bilateral  hull  symmetry,  and  the  waterplane,  in  accordance 
with  the  zero  Froude  number  assumption.  With  this  symmetry,  only  the  quadrant  in  which 
y'  and  z'  are  both  positive  is  considered  directly  in  the  potential  flow  calculations,  the 
other  three  quadrants  being  accounted  for  as  'images'  in  the  symmetry  planes. 


3  CALCULATION  OF  THE  DISCRETE  CARTESIAN  POTENTIAL  FLOW 


The  calculation  of  the  incompressible  potential  flow  velocities  for  flow  around  a  ship 
hull  consists  of  solving  Laplace's  equation  subject  to  given  boundary  conditions.  It  has 
long  been  known  that  this  can  be  done  efficiently  using  boundary  integral  methods, 
effectively  reducing  the  three-dimensional  problem  to  a  two-dimensional  one.  In  order  to 
evaluate  the  integrals  over  arbitrary  boundaries,  recourse  had  to  be  made  to 
discretisation,  which  lead  to  the  historical  development  of  panel  methods  exemplified  by 

the  early  work  of  Hess  and  Smith5. 

The  method  used  in  present  work  is  essentially  as  formulated  by  Hess  and  Smith5. 
For  a  rigorous  development  of  the  theory,  the  reader  is  referred  to  the  article  by  Hunt6  on 
boundary  integral  methods.  The  principal  steps  are  outlined  below. 

Figure  4  represents  a  three-dimensional  domain  ft  with  boundary  S  on  which  the 
normal  n  is  directed  into  ft.  P  is  any  point  in  the  domain,  or  flow  region.  The  unknown 
velocity  potential  in  the  domain,  4>,  satisfies  Laplace's  equation  in  this  region: 

V2#  =  0  in  ft  (3.1) 

It  is  convenient  to  split  $  into  two  terms: 

$  =  $s  +  $«,  (3.2) 

$s  is  an  inner  solution  which  induces  velocities  which  vanish  at  infinity.  Since  the  flow 
has  no  circulation  and  there  is  no  free  surface,  an  appropriate  form  for  $s  is  that  derived 

from  a  simple  (Rankine)  source7  distribution  on  the  inner  boundaries;  the  boundary 
condition  that  the  the  velocities  induced  by  $ s  vanish  at  infinity  is  then  automatically 
satisfied. 


*s 


1 

47TT 


(3.3) 


47TT  is  the  solid  angle  subtended  by  the  boundary  at  the  point  at  which  the  potential  is 
evaluated.  If  this  point  is  within  the  domain  ft,  the  solid  angle  will  be  47T;  if  the  point  is  on 
S,  the  solid  angle  will  be  less  than  47T. 

accounts  for  the  boundary  condition  that  the  velocity  be  uniform  at  infinity: 

$«,  *  -x*Vw  (3.4) 

where  V.  is  the  free  stream  velocity.  Notice  that  both  $s  and  satisfy  Laplace's 
equation  separately.  The  boundary  condition  satisfied  at  S  is 

n-V$  *  0  (3.5) 


which  may  be  expressed  as 
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n-V$s  =  n’V,,,  on  S  (3.6) 

The  problem  has  now  been  reduced  to  determining  the  unknown  source  density  O' 
such  that  the  boundary  condition  of  equation  (3.6)  is  satisfied.  Since  the  surface  integral 
in  equation  (3.3)  can  only  be  evaluated  analytically  for  certain  regular  shapes,  in  general 
the  problem  must  be  discretized  and  the  surface  integral  replaced  by  a  summation  as 
illustrated  by  Figure  5.  The  surface  is  discretized  into  a  number  of  panels  each  with  an 
associated  control  point  at  which  the  boundary  condition  is  enforced.  The  source 
distribution  is  uniform  over  the  extent  of  each  panel.  The  inner  potential  $s  may  then  be 
expressed  as 


•.-!>*  «-7> 

j 

where 


where  4>Sj  is  the  contribution  to  the  potential  from  a  unit  source  distribution  on  the 
panel.  The  velocity  at  the  1th  control  point  X|  in  the  domain  is 


V, 


(3.9) 


where 

•u  ■  -V* j(x,).  (3.10) 

a,j  is  the  velocity  induced  at  the  ith  control  point  due  to  a  unit  source  distribution  on  the 
jm  panel.  They  are  called  the  influence  coefficients. 

The  boundary  condition  of  equation  (3.6)  evaluated  at  the  Ith  control  point  is 

nrv» = 2  n,'a,j<rj  (3,11) 

j 

The  oys  are  evaluated  by  inverting  equation  (3.11).  The  tangential  velocity  components 
and  velocities  in  the  field  are  then  obtained  directly  from  equation  (3.9).  Pressures  may 
be  calculated  with  Bernoulli's  equation7. 

In  the  discrete  model,  each  influence  coefficient  contains  the  surface  Integral  of 
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equation  (3.3)  evaluated  over  the  j^1  panel  only.  For  the  present  work,  the  panel  is 
assumed  to  be  planar,  and  the  source  distribution  to  be  constant  on  its  surface.  This  is 
the  first-order  approximation  and  results  in  relatively  simple  analytic  expressions  for  the 
surface  integral. 


4  COMPRESSIBILITY  CORRECTION  FOR  POTENTIAL  FLOW 

Where  the  calculation  is  being  done  for  comparison  with  data  from  a  model  in  a  wind 
tunnel,  compressibility  of  the  fluid  must  be  taken  into  account  when  the  local  Mach 
number  becomes  too  large.  A  rule  of  thumb  is  that  corrections  should  be  made  when  the 
free-stream  Mach  number,  Mw,  exceeds  0.3.  Such  compressible  flows  are  governed  by 

the  Prandtl-Glauert  equation, 


(1  -  m£ 


(4.1) 


so  long  as  the  local  Mach  number  does  not  exceed  1 .0.  The  Prandtl-Glauert  equation  can 

be  reduced  to  Laplace's  equation  by  the  application  of  Gothert's  rule8,  a  coordinate 
transformation  which  stretches  the  hull  in  the  x1  direction.  The  calculated  incompressible 
velocity  potential  is  then  transformed  back  to  compressible  coordinates. 


5  SPLINING  THE  DISCRETE  CARTESIAN  POTENTIAL  FLOW 

In  this  section  the  technique  used  by  the  program  VELSPL  to  spline  the  discrete 
cartesian  potential  flow  will  be  described. 


5.1  REPRESENTATION  OF  THE  VELOCITY  VECTOR  IN  HULL  COORDINATES 

Since  the  hull  coordinates  are,  in  general,  non-orthogonal,  there  are  several  ways  in 
which  the  flow  vectors  can  be  represented  in  these  coordinates.  The  most  convenient  is 
to  provide  the  covariant  or  the  contravariant  components  of  the  vector.  Although  the 
covariant  components  of  the  vectors  are  the  most  straightforward  to  calculate  the 
representation  in  terms  of  contravariant  components  has  been  chosen  for  two  reasons. 

1)  It  is  required  that  the  keel  and  the  waterline  are  streamlines.  These 
conditions  are  most  easily  enforced  using  the  contravariant  representation 
(see  Section  5.3). 

2)  The  calculation  of  streamlines  is  simplified  by  the  use  of  the  contravariant 
representation  (see  Appendix  C). 

The  contravariant  components  of  a  vector  in  the  non-orthogonal  surface  coordinate 
system  can  be  determined  If  the  relationship  between  these  coordinates  and  the 
cartesian  coordinate  system  is  known.  The  calculation  method  is  outlined  below.  Define 
two  linearly  independent  hull-surface  vectors,  ek,  k  *  1,2,  by 


.*•  -•> .' 


.>  .■>  > , 
v.v.v\,w.\ 
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where,  for  convenience,  the  position  coordinate  (x,s)  has  been  denoted  (xvx2).  Two 
additional  linearly  independent  hull-surface  vectors,  ek,  may  then  be  defined  such  that 

ek-en  =  $kn  (5.2) 


where  6kn  is  the  Kronecker  delta.  Simple  algebraic  manipulation  then  gives 

(er>'®k)*k  *  ®n  (5-3> 

where  summation  over  repeated  indices  Is  implicit.  The  ek  may  therefore  be  obtained  from 
the  ek  by  the  inversion  of  the  2x2  matrix  e^e*.  The  contravariant  components,  Vk  of  a 
vector,  V,  are  defined  by 

Vk  =  v-ek  (5.4) 


Since  the  derivatives  dX/dx  may  be  calculated  easily  using  the  spline  coefficients 
supplied  by  SMHULL1  the  calculation  of  the  contravariant  components  at  the  control  points 
is  straightforward. 


5.2  SPLINING  THE  CONTRAVARIANT  VELOCITY  COMPONENTS 

Once  the  contravariant  components  of  the  velocity  vectors  are  known  at  the  control 
point'-,  these  values  are  interpolated  using  a  smoothing  spline.  The  resulting  two 
dimensional  spline  may  be  used  to  determine  the  value  of  the  contravariant  velocity 
components  at  any  point  on  the  hull.  The  method  of  spiining  is  very  similar  to  that  used 

for  spiining  the  cartesian  coordinates  In  the  program  SMHULL1,  and  so  will  be  treated  in 
less  detail  here. 

The  use  of  smoothing  spline  rather  than  splines  which  interpolate  the  data  is 
prompted  by  two  observations:  that  splines  can  be  very  sensitive  to  the  placement  of  the 
data  points  and  the  overall  spline  fit  can  often  be  improved  dramatically  by  allowing  the 
spline  to  vary  very  slightly  from  the  data  points;  and  that  the  velocities  calculated  by 
EN967  cannot  be  considered  without  error,  so  that  improvements  of  the  splines  within  the 
tolerances  of  the  errors  expected  from  EN967  can  only  improve  resulting  calculations.  It 
should  also  be  emphasized  that  the  user  has  control  over  the  amount  of  smoothing  done 
by  the  spline  (see  Appendix  B.4.1). 
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5.2.1  REPRESENTATION  OF  THE  VELOCITY  IN  TERMS  OF  B-SPLINES 

The  contravariant  components  of  the  velocity  V1  and  V2  are  represented  as  two- 
dimensional  splines  in  terms  of  the  independent  variables  x  and  s.  The  control  points 
must  form  a  rectangular  grid  in  the  non-orthogonal  coordinates,  (x,s);  that  is,  the  control 
points  must  be  of  the  form  {(Xp,sq):  p=1  ,NX;  q=1,Ny).  This  is  ensured  by  the  program 
EN967S  when  it  chooses  the  control  points. 


Nx  N, 

Vk(X,s)  =  £E*«  BP.kx(x)Bq.kI(s) 

p=1  q=1 


(5.5) 


where  Bp  k  (x),  p=1,Nx,  is  a  B-spline  of  order  kx  corresponding  to  a  knot  sequence 

t£x\  p=1  ,Nx+kx,  and  similarly  for  Bqt^(s).  Notice  that  the  orders  of  the  B-splines  in  the  x 

and  s  directions  need  not  be  the  same,  nor  need  these  be  the  same  as  the  orders  of  the 
knot  sequences  used  in  the  SMHULL  spline  representation,  despite  the  similarity  in 
notation  used  here  and  in  Section  2. 

For  the  purposes  of  computation,  it  is  most  convenient  to  express  equation  (5.5)  in 
the  following  form: 

N. 

V'^x.s)  «  ^  ukq(x)BPfk>(s)  ,  (6.6) 

q«1 


* 

Mkq(x)  «  2  XlWWX)- 


5.2.2  EVALUATION  OF  THE  SPLINE  COEFFICIENTS 

The  evaluation  of  the  spline  coefficients  can  be  broken  into  two  steps: 

1 )  On  each  curve  of  constant  x,  the  contravariant  velocity  components  are  splined 
in  the  s  direction  using  the  knots  q*1,N,+k,.  This  yields  the  spline 
coefficients  ykq(Xp),  where  Xp  is  the  p-th  curve  of  constant  x. 

2)  The  spline  coefficients  ukq(Xp)  are  themselves  splined  in  the  x  direction  using 

the  knots  t£x\  p«1,Nx+kx.  This  yields  the  required  spline  coefficients  ^kw.  For 
this  step  to  be  possible  it  is  necessary  that  the  same  number  of  B-splines  be 
used  to  spline  each  offset  station  and  that  the  knots  vary  smoothly  as  functions 
of  x.  Here,  the  same  knots  are  taken  for  every  offset  station,  reducing  the 
necessary  storage  space. 
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The  knots  for  the  B-splines  may  be  input  directly  by  the  user  or  default  values  may 
be  accepted  (see  Appendix  ). 

Techniques  for  constructing  smooth  spline  curves  approximating  given  data  have 
been  described  in  detail  by  de  Boor3  and  by  Hally1,4.  The  subroutine  BSMTH  described  by 

Hally4  has  been  chosen  for  use  in  VELSPL.  It  determines  a  smooth,  accurate 
approximation  to  the  given  data,  with  good  control  over  the  'wiggles'  which  plague  other 
spline  methods. 


5.3  END  POINT  CONDITIONS  FOR  THE  SPLINES 

The  velocities  calculated  by  EN967  are  given  at  the  control  points  of  the  planar 
panels.  Since  the  control  points  are  near  the  centres  of  the  panels,  there  are  no 
velocities  given  at  the  boundaries  of  the  the  hull,  keel  and  the  waterline.  However,  the 
symmetry  of  the  hull  requires  that  the  keel  and  the  waterline  be  streamlines,  so  the 
splines  of  the  velocity  data  must  take  this  Into  account.  Since  the  keel  and  the  waterline 
are  curves  of  constant  hull  coordinate  s,  the  end  point  condition  applied  to  the  splines 

may  be  restated  as  follows.  Since  et  is  a  vector  parallel  to  the  curve  of  constant  s,  e2  is 
a  vector  perpendicular  to  the  curve  of  constant  s  (see  equation  5.2).  Hence,  the 
condition  that  a  curve  of  constant  s  be  a  streamline  is 

V-e2  =  V2  =  0  (5.8) 

Hence,  the  streamline  condition  to  be  applied  to  the  splines  reduces  to  a  constraint  on 
the  value  of  V2  at  the  boundary.  This  constraint  is  implemented  most  easily  by  requiring 
that  the  first  k,  knots  in  the  s-direction,  p=1,ks  lie  at  s  =  0  (the  waterline)  and  that 
the  last  ks  knots  p=N4+1,Ns+k#,  lie  at  s  =  1  (the  keel).  One  can  then  show3  that 

liVx)  *  V2(x,0)  ;  m\is(x)  =  V2(x,1).  (5.9) 


Thus,  the  waterline  and  keel  can  be  constrained  to  be  streamlines  simply  by  requiring  that 
X21q  »  0  and  that  X2Nj()  =  0. 


6  CONCLUDING  REMARKS 

Incompressible  potential  flow  gives  an  adequate  description  of  the  flow  around  a 
ship  hull  for  many  purposes  including  the  computation  of  surface  velocities  for  a  boundary 
layer  calculation.  Details  of  a  numerical  method  for  calculating  this  flow  have  been 
presented  together  with  a  method  for  interpolating  the  output  data  to  make  it  more 
suitable  for  use  in  further  hydrodynamic  computations. 

The  appendices  to  this  memorandum  provide  user  guides  to  computer  programs 
which  calculate  potential  flow  and  which  spline  the  resulting  surface  velocities.  These 
programs  form  part  of  a  total  package  for  calculating  the  development  of  ship  hull 

boundary  layers,  the  other  parts  of  which  are  documented  by  Hally1,2. 
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HULL  DATA 


DIGITIZE  HULL  DATA 

GENERATE  HULL  B-SPLINE 
COEFFICIENTS 

SET  UP  DATA  FOR  EN967 


CALCULATE  POTENTIAL 
FLOW  VELOCITIES 


GENERATE  VELOCITY 
B-SPLINE  COEFFICIENTS 


’CALCULATE  BOUNDARY 
LAYER 


#1:  DESCRIBED  IN  REFERENCE  1 
#:  DESCRIBED  IN  THIS  MEMORANDUM 
#2:  DESCRIBED  IN  REFERENCE  2 


Figure  1 .  Flow  Chart  of  Boundary  Layer  Calculation 
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Figure  5.  Equivalence  of  Potential  Flow  and  Discrete  Model 
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Appendix  A 

USER'S  GUIDES  FOR  THE  EN967  PROGRAMS 


The  program  suite  EN967  was  obtained  from  the  DeHavifiand  Aircraft  Company  of 
Canada  Ltd.  It  models  lifting  surfaces  and  therefore  includes  line  vortex  singularities  in 
addition  to  sources.  The  original  'front  end’  was  optimized  for  the  representation  of 
aircraft  geometries  and  consisted  of  two  programs: 

EN967A  -  sets  up  total  geometry  from  an  input  description  of  cross-sections 
etc.,  in  terms  of  non-pianar  panels  (ie.:  panels  whose  corner  points 
lie  on  the  true  surface). 

EN967B  -  calculates  equivalent  planar  panels  from  which  influence 
coefficients  can  be  calculated.  Locates  control  points. 

For  lifting  surfaces,  EN967A  and  EN967B  also  perform  analogous  functions  for  the 
vortices. 

For  ship  hulls  the  above  two  programs  are  replaced  by  a  single  program,  EN967S, 
which  takes  as  input  data  the  hull  B-spline  coefficients  on  a  data  file  created  by 

SMHULL1.  EN967S  mimics  the  two  programs,  producing  disc  files  and  lineprinter  output  in 
the  original  format,  the  latter  being  labelled  EN967A  or  EN967B  equivalent.  In  addition, 
EN967S  writes  the  input  files  EN967C.DAT,  EN967D.DAT  and  EN967E.DAT  required  to 
complete  the  potential  flow  calculation,  so  that  the  remaining  steps  are  virtually 
transparent  to  the  user. 

The  EN967  programs  save  significant  intermediate  results  on  binary  disc  files  called 
NAME*.DAT,  where  *  is  an  integer  between  7  and  30.  A  list  of  these  files  and  their 
functions  is  given  in  Table  A.1.  Note  that  although  no  vortices  are  used  in  the  ship  hull 
calculations,  EN967S  creates  the  file  NAME7.DAT  containing  the  vortex  data 
nevertheless,  so  that  EN967C  will  not  fail  because  it  can't  find  the  input  file.  NAME7.DAT 
is  empty  if  created  by  EN967S.  Similarly,  EN967C  creates  empty  files  NAME15.DAT  and 
NAME17.DAT  for  use  in  EN967D  and  EN967E  respectively. 

The  NAME*  files  generate  a  considerable  storage  requirement  if  a  large  number  of 
hull  panels  is  used.  Storage  on  the  DREA  DEC  2060  is  measured  in  pages;  a  page 
consists  of  512  16-bit  words.  Table  A.1  gives  the  storage  requirement  for  a  1000  panel 
example.  Note  that,  as  a  rule  of  thumb,  NAME*  file  storage  is  proportional  to  the  square 
of  the  number  of  panels. 


A.1  USER'S  GUIDE  FOR  EN967S 

Program  EN967S  generates  the  non-planar  and  equivalent  planar  panel  geometry 
from  B-spline  coefficients  supplied  by  program  SMHULL.  The  program  initially  distributes 
non-planar  panels  so  that  their  corner  points  are  equally  spaced  in  (x,s)  coordinates.  The 
number  of  panels  in  each  direction  is  requested  by  the  user.  The  user  may  then  modify 
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the  panel  distribution  so  that  it  is  more  optimal  for  the  potential  flow  calculation:  the 
lengthwise  (x)  distribution  is  altered  so  that  regions  with  great  curvature  have  a  higher 
panel  density  than  regions  with  a  lesser  curvature.  This  is  achieved  by  making  panel 
lengths  inversely  proportional  to  the  second  x-derivative  of  the  arclength  in  the  s 
direction.  Adjacent  panel  lengths  are  constrained  so  that  their  ratio  does  not  exceed 
1.5.  This  modified  panel  distribution  will  have  more  panels  in  the  x  direction  than  were 
initially  requested. 

Figure  A.1  illustrates  the  principal  geometrical  features  of  a  non-planar  panel 
(with  curvatures  greatly  exaggerated)  and  the  equivalent  planar  panel.  The  non-planar 
panel  corner  points,  P,,  i=1,4,  and  the  panel  control  point,  Pcp,  are  determined  from  the 
spline  coefficients  supplied  by  SMHULL,  as  are  the  gradients  of  the  surface  coordinates 
with  respect  to  the  cartesian  coordinates,  x  and  s.  The  panel  normal,  n,  is  directed  along 
the  vector  product  of  these  gradients.  The  equivalent  planar  panel  is  located  by  point  PN 
on  the  normal,  where  PN  is  the  average  of  the  control  point  and  the  non-planar  comer 
points.  Corner  points  of  the  planar  panel,  P',,  i*1,4,  are  projections  o'  the  non-planar 
panel  corner  points. 

The  panel  control  point  must  be  within  the  flow  domain,  or  on  its  boundaries.  For  a 
concave  non-planar  panel  the  locations  of  Pqp  and  PN  are  reversed  from  the  configuration 
shown  in  Figure  A.1,  placing  the  control  point  outside  the  domain.  In  this  case,  PN  is  co¬ 
located  with  Pcp. 

EN967S  calculates  a  number  of  additional  properties  of  the  planar  panel:  side 
lengths  and  slopes,  and  zeroth  (area),  first  and  second  moments  of  inertia.  These  are 
required  to  determine  the  influence  coefficients  in  EN967C. 


A.  1.1  OPERATION  OF  EN967S 

To  run  EN967S  on  the  DREA  DEC-20  computer 

EXECUTE  PS: <I1ACKAY>EN367S,  PS: <HALLY>SMSUBS/l_IB,  PS: <HALLY>GETURD/LIB, 

BSPL I N : HLL  YSP/L I B ,  BSPL I N : BSPLN7 /L I B 

Program  SMHULL  supplies  an  ASCII  free-format  file  (except  for  record  1)  which  is 
read  by  EN967S.  The  contents  of  the  file  are  outlined  in  Table  A.2.  In  addition,  there  is  a 
small  amount  of  dialog  with  the  user  at  execution  time:  this  is  described  below.  Items 
labelled  'prompt1  and  'comment1  are  generated  by  the  program;  Items  labelled  'response' 
are  entered  by  the  user. 

Prompt  1:  ENTER  P  FOR  LPT:  OUTPUT  -  OTHERWISE  ON  EN9S7S.0UT 

Response:  The  character  'P'  or  'p'  directs  the  program  output  tables  to  the  printer. 

Any  other  response  (including  <CR>)  causes  the  output  to  be  saved  on 
file  EN967S.0UT. 


Prompt  2:  Enter  name  of  file  containing  hull  spline  coefficients  »> 
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Response:  FILDUM:  the  name  (extension  optional,  the  default  is  .DAT)  of  the  SMHULL 
data  file.  A  caret  (a)  aborts  the  run. 

Comment:  Reading  data  from  FILDUM  ...  Done 
...HULL  DATA  FILE  TITLE  IS  : 

TYTLE  (up  to  60  characters) 

Prompt  3:  ENTER  NO.  OF  PANELS  IN  X  DIRECTION  ... 

...  IE:  NO.  OF  SECTIONS  BOU  TO  STERN  (MAX  43) 

Response:  NXC:  a  positive  integer,  less  than  or  equal  to  limit  NXCM  (currently  49). 

Prompt  4:  ENTER  NO.  OF  PANELS  AROUND  GIRTH  (MAX  49) 

Response:  NSC:  a  positive  integer,  less  than  or  equal  to  limit  NSCM  (currently  49). 

Prompt  5:  ENTER  G  TO  USE  VARIABLE  GRID  OPTION 

Response:  The  character  'G'  or  'g'  results  in  a  grid  adjustment  as  described  in  the 
general  discussion  of  the  program.  The  number  of  panels  in  the  x- 
direction  will  be  increased.  Any  other  response  (including  <CR>)  will 
cause  this  option  to  be  ignored. 

Comment  :  NO.  OF  X-PANELS  IS  NXC  (where  NXC  is  revised  as  above) 

Prompt  6:  ENTER  T  FOR  GOTHERT  TRANSFORMATION 

Response:  The  character  *T'  or 't'  indicates  that  compressibility  corrections  will  be 
required  and  Prompt  7  will  be  given.  Any  other  response  (including 
<CR»  will  result  in  the  Mach  number  being  set  to  zero. 

Prompt  7:  ...  MACH  NUMBER  ? 

Response:  AMACH:  freestream  Mach  number. 

In  general,  an  unacceptable  response  will  result  in  a  self-explanatory  error 
message  and  a  repeat  of  the  appropriate  prompt.  Note  in  particular  that  there  is  a  limit 
on  the  number  of  panels  which  can  be  accommodated  by  EN967.  This  limit  is  1200  in  the 
standard  version  of  the  program.  If  the  limit  is  exceeded,  either  by  choice  of  NXC  and 
NSC,  or  inadvertently  when  using  the  grid  adjustment  option,  an  error  message  is  given 
and  the  program  returns  to  Prompt  3. 

Execution  of  EN967S  is  quite  fast:  25  seconds  of  CPU  time  for  1000  panels 
(25x40)  without  grid  adjustment,  and  35  seconds  with  grid  adjustment  (17x40  ■* 
25x40).  CPU  time  is  approximately  proportional  to  the  number  of  panels. 


A.  1.2  OUTPUT  FROM  EN967S 

EN967S  produces  tables  of  results  and  disc  files  for  the  use  of  subsequent  program 
units.  The  tables  are  either  printed  immediately  or  stored  on  a  file  EN967S.OUT, 
according  to  the  response  to  Prompt  1  in  the  previous  section.  Disc  file  output  consists 
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of  ASCII  formatted  files  EN967C/D/E.DAT  and  the  binary  files 
NAME7/9/1 2/1  7/1 8/24.DAT. 

The  results  tables  from  EN967S  include 

-  date  and  time  of  run 

-  title,  with  date  and  time  appended;  this  version  of  the  title  is  passed  on  to 
subsequent  EN967  units  and  provides  a  unique  run  identifier 

-  NXC.NSC:  numbers  of  panels  in  x  and  s  direction 

-  length,  half-breadth  and  depth  of  hull 

-  x'.y'.z1  coordinates  of  P1.P2.P3.P4.  the  non-planar  panel  corner  points 
("EN967A  equivalent  output") 

-  x',y',z'  elements  of  the  vectors  T.,  and  T2  defining  the  planar  panel  and  the 
normal  vector  n,  and  x'.y'.z*  coordinates  of  the  control  point  Pcp  and  the 
planar  panel  definition  point  PN  (refer  to  Figure  A.1  for  definitions). 
Vectors  T-,  and  T2  join  the  corners  of  the  planar  panel  and  are  normal  to  n 
("EN967B  equivalent  output"). 

The  content  and  format  of  output  files  EN967C/D/E.DAT  are  described  in  the  input 
sections  of  the  appropriate  program.  The  purpose  of  files  NAME7/9/ 12/1  7/1 8/24.DAT  is 
summarized  in  Table  A.1. 

The  user  who  does  not  wish  to  be  concerned  with  the  details  of  programs 
EN967C/D/E  may  proceed  directly  to  Section  A.5  for  the  commands  to  run  these  units. 


A.2  USER'S  GUIDE  FOR  EN967C 

Program  EN967C  calculates  the  elements  atj  of  the  influence  coefficient  matrix.  It 
remains  essentially  as  coded  at  DeHavilland,  except  that  the  self-influence  calculation  is 
modified  so  that  the  control  point  need  not  lie  on  the  planar  panel. 

The  influence  of  a  planar  source  panel  is  calculated  with  the  formulae  given  by  Hess 

and  Smith5;  these  formulae  will  not  be  repeated  here.  A  characteristic  dimension  of  the 
panel,  the  longest  diagonal,  is  chosen,  and  the  ratio  of  the  distance  to  the  panel  and  the 
characteristic  dimension  is  determined.  For  computational  efficiency,  the  influence 
coefficient  is  then  given  by 

-  far-influence,  if  the  ratio  is  greater  than  6.  A  point  source  approximation 
is  used. 

-  mid-influence,  if  the  ratio  is  between  2.45  and  6.  A  multipole  expansion  of 
the  source  density  is  used. 
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*  near-influence,  if  the  ratio  is  less  than  2.45.  The  surface  integral  over  the 
planar  panel  is  evaluated  exactly. 

To  avoid  numerical  difficulties  with  self-influence  when  Pcp  and  PN  are  not  co-located,  the 
near-influence  formula  for  the  normal  component  is  modified  by  an  expansion  with  respect 
to  the  distance  of  Pcp  above  the  panel.  If  Pcp  and  PN  are  co-located,  the  normal 
evaluation  of  self-influence  is  used. 

As  only  one  quadrant  (both  y1  and  z‘  positive)  of  the  double  hull  geometry  is 
explicitly  accounted  for  (Figure  3),  each  influence  coefficient  represents  the  sum  of 
terms  from  all  four  quadrants. 


A.2.1  OPERATION  OF  EN967C 

To  run  EN967C  on  the  DREA  DEC-20  computer: 

EXECUTE  PS: <MACKAY>ENS67C 

EN967C  requires  an  ASCII  formatted  data  file  EN967C.DAT,  and  binary  files 
NAME7 /9/ 24.DAT,  created  in  EN967S.  The  contents  of  EN967C.DAT  are  outlined  in  Table 
A.3.  The  default  values  shown  are  supplied  by  EN967S. 

A  case  with  1 000  panels  and  two  planes  of  symmetry  takes  about  34  CPU  minutes; 
execution  time  is  approximately  proportional  to  the  number  of  panels  squared. 


A.2.2  OUTPUT  FROM  EN967C 

EN967C  prints  one  page  of  output: 

-  date  and  time  of  run 

-  KODE,  KXOUT  (see  Table  A.3) 

-  title 

-  counters  which  indicate  how  many  times  each  method  for  calculating 
source  influence  is  used:  KEX  for  self-  and  near-influence,  KflP  for  mid¬ 
influence  and  KPS  for  far-influence.  The  total  of  these  is  four  times  the 
square  of  the  number  of  panels  for  two  planes  of  symmetry 

-  message:  "PART  C  COMPLETED" 


In  addition,  if  KXOUT  -  1,  the  matrix  of  influence  coefficients  is  printed  in  1  PEI  2.4  format. 
The  function  of  binary  disc  files  NAME  10/1 1/1 3/1 5.DAT  generated  by  EN967C  is 
summarized  in  Table  A.1. 
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A.3  USER’S  GUIDE  FOR  EN967D 

Program  EN967D  solves  the  set  of  simultaneous  equations  represented  by  equation 
(3.9).  It  remains  essentially  as  coded  at  DeHavllland. 

The  iterative  method  described  by  Labrujere  et  al.s  is  used.  The  influence  matrix  is 
'blocked1,  or  subdivided.  Since  matrix  entries  near  the  diagonal  are  dominant  due  to  large 
self-  and  near>influence  terms,  solving  the  submatrices  on  the  diagonal  alone  gives  a  fair 
approximation  to  the  inverse  of  the  full  matrix.  A  correction  term  is  formulated  and  the 
calculation  is  iterated  until  the  deviation  from  the  required  tangential  flow  boundary 
condition  (the  residual)  is  acceptable.  Less  than  10  iterations  is  usually  sufficient  for 
convergence.  The  calculation  is  terminated  after  30  iterations  if  convergence  has  not 
occurred;  this  condition  usually  indicates  that  mixed  NAME*.DAT  files  from  different  runs 
are  being  used  since  the  influence  matrix  is  invariably  well  conditioned  for  ship  hulls.  A 
Gauss-Seidel  algorithm  is  used  to  invert  the  submatrices. 


A.3.1  OPERATION  OF  EN967D 

To  run  EN967D  on  the  DREA  DEC-20  computer 
EXECUTE  PS: <MACKAY>EN9S7D 

EN967D  requires  an  ASCII  formatted  data  file  EN967D.DAT,  and  binary  files 
NAME1 0/1 5.DAT,  created  in  EN967C.  The  contents  of  EN967D.DAT  are  outlined  in  Table 
A.4.  The  defaults  are  set  by  EN967S. 

A  case  with  1000  panels  takes  about  4.5  CPU  minutes;  execution  time  is 
approximately  proportional  to  the  number  of  panels. 


A.3.2  OUTPUT  FROM  EN967D 

EN967D  produces  the  following  lineprinter  output  according  to  the  default  value  of 
NSB  (see  Table  A.4); 

-  matrix  blocking  information 

-  a  number  of  items  not  relevant  to  use  with  EN967S,  including  vortex,  angle 
of  attack  and  submatrix  inversion  data  (these  can  be  significant  for  lifting 
surface  applications) 

-  sum  of  squares  of  residuals  for  each  submatrix 

The  angles  of  attack  (alpha)  and  sideslip  (beta)  referred  to  in  this,  and  EN967E,  output 
are  constrained  to  be  zero  with  EN96 7S  since  the  freestream  flow  must  be  tangential  to 
both  planes  of  symmetry  (Figure  3). 

If  NSB  is  given  the  alternative  value  noted  in  Table  A.4,  a  detailed  table  is  printed 
including  source  densities  and  residuals  for  the  last  iteration. 


*«■  •  *•  <*  .*  •  •  •  ■»  .  _  .  *.  « 

•  •  *.*  .•  .  •*.  • 


.  •  •  ’ -  *  .  *  •  . 


■ v>*. ^ 
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The  function  of  binary  disc  file  NAME16.DAT  generated  in  EN967D  is  summarized  in 
Table  A.1. 


A.4  USER'S  GUIDE  FOR  EN967E 

Program  EN967E  calculates  surface  velocities  using  equation  (3.9).  Primarily  of 
interest  in  aerodynamic  applications  rather  than  the  task  addressed  here,  are  a  number  of 
derived  quantities: 


pressure  coefficients, 


-  force  and  moment  coefficients 


The  potential  flow  about  a  symmetrical  nonlifting  body  gives,  within  numerical  error,  zero 
total  forces  and  moments. 


A.4.1  OPERATION  OF  EN967E 

To  run  EN967E  on  the  DREA  DEC-20  computer 
EXECUTE  PS: <HACKA Y >EN967E 

EN967E  requires  an  ASCII  formatted  data  file  EN967E.DAT,  and  binary  files 
NAME  11/12/13/16/ 24. DAT,  created  in  EN967S/C/D.  The  contents  of  EN967E.DAT  are 
outlined  in  Table  A.5.  The  defaults  are  set  by  EN967S. 

A  case  with  1 000  panels  and  two  planes  of  symmetry  takes  about  9.5  CPU  minutes; 
execution  time  is  approximately  proportional  to  the  square  of  the  number  of  panels. 

If  the  default  1 1  LPT  (see  Table  A.5)  has  been  accepted  and  tabular  output  is 
required, 

PRINT  EEOUT/FILE: FORTRAN 


A.4.2  OUTPUT  FROM  EN967E 

EN967E  saves  on  EE0UT.DAT,  or  prints,  the  following  lineprinter  output  according  to 
the  default  value  of  1 1  LPT : 

-  run  date  and  time 

-  title  and  Mach  number 

-  for  each  control  point:  x',y',z'  coordinates,  x'.y'.z1  components  of  flow 
velocity,  source  density  and  pressure  coefficient 

-  force  and  moment  coefficients  by  block  ('section')  and  totals 

If  EE0UT.DAT  is  saved,  a  single  page  of  lineprinter  output  is  generated  indicating 
completion  and  reminding  the  user  to  print  EEOUT.DAT  (if  required). 


*  '»*.*"  .* '*'■_■**  i**  i*' .*■  ,  *  .**  A  A  *»'  >'  *.'  »B*  ^ 


.  V 
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EN967E  saves  velocities  and  pressure  coefficients  on  NAME30.DAT  for  the  use  of 
VELSPL,  by  default. 


A.S  SUMMARY  OF  OPERATING  PROCEDURE  FOR  ENS67C/D/E 

In  the  preceeding  sections,  programs  EN967C,  EN967D  and  EN967E  have  been 
outlined.  However,  since  EN967S  prepares  the  necessary  input  files,  the  casual  user 
need  not  bother  with  these  details,  but  simply  invoke,  on  the  DREA  DEC-20  system, 

EXECUTE  PS: <MACKAY>EN967C 
EXECUTE  PS: <MACKAY>EN967D 
EXECUTE  PS: <MACK A Y >EN967E 
PRINT  EEOUT/F I LE: FORTRAN 

Note  however,  that  user  directories  will  generally  be  too  small  to  accommodate  the 
NAME*. DAT  files,  so  that  these  units  will  have  to  be  run  on  DREA:<DREA>.  The  NAME*.DAT 
files  will  remain  in  the  host  directory  until  deleted. 

Disc  storage  requirements  are  summarized  in  Table  A.1  and  execution  times  in 
Figure  A.2. 
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Table  A.1 

TEMPORARY  DISC  FILE  STORAGE  FOR  EN9S7 


FILENAME 

CREATED 

IN 

USED 

IN 

CONTAINS 

STORAGE  IN 
PAGES,  FOR 
1 000  PANELS 

NAME7.DAT 

EN967S 

EN967C 

vortex  data 

0 

NAME9.DAT 

EN967S 

EN967C 

planar  panel  data 

79 

NAME10.DAT 

EN967C 

EN967D 

normal  induced  velocities 

1975 

NAME  1 1  .DAT 

EN967C 

EN967E 

source  influence  matrix 

5917 

NAME12.DAT 

EN967S 

EN967E 

control  point  data 

18 

NAME13.DAT 

EN967C 

EN967E 

normal  vortex  induced  velocities 

0 

NAME  1 5.DAT 

EN967C 

EN967D 

vortex  influence  matrix 

0 

NAME16.DAT 

EN967D 

EN967E 

singularity  strengths 

2 

NAME17.DAT 

EN967S 

VELSPL 

surface  derivatives 

8 

NAME  1 8.DAT 

EN967S 

non-planar  panel  data 

28 

NAME24.DAT 

EN967S 

EN967C/E 

geometry  counters 

3 

NAME30.DAT 

EN967E 

VELSPL 

velocities  and  pressure  coefs. 

18 

NOTES 

1)  Files  associated  with  vortex  data  are  empty  (NAME7/13/15). 

2)  The  storage  for  NAME  10  and  NAME  11  is  proportional  to  the  square  of  the  number  of 
panels;  for  other  files  it  is  proportional  to  the  number  of  panels. 

3)  NAME  1 8  is  included  to  make  EN967S  output  equivalent  to  that  of  EN967A;  the  file  is 
not  used  for  boundary  layer  calculations.  New  graphics  software  is  being  planned 
which  will  use  this  file. 
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Table  A.2 

CONTENTS  OF  EN967S  INPUT  DATA  FILE 


Record 

Variable 

Format 

Purpose 

1 

TYTLE 

A60 

60  character  title 

UNITS 

A4 

units  used 

2 

SLNGTH 

free 

LBP 

BREDTH 

free 

maximum  haif-breadth 

DEPTH 

free 

maximum  draft 

3 

NST 

free 

number  of  stations 

XS(IX) 

free 

station  numbers,  IX-1,NST 

4 

KS 

free 

order  of  B-splines  in  s-direction 

NS 

free 

number  of  B-splines  in  s-direction 

TSU) 

free 

s-knot  locations,  I  -1 ,  KS+NS 

5 

KX 

free 

order  of  B-splines  in  x-direction 

NX 

free 

number  of  B-splines  in  x-direction 

TX(I) 

free 

x-knot  locations,  I  -1 ,  KX+NX 

6 

BCOEFY (IS, IX) 

» 

B-spline  coefficients  a(j  and  /3n 

BCOEFZ (IS, IX) 

free  IS-l.NS;  IX-l.NX 

Table  A.3 

CONTENTS  OF  EN9G7C.DAT 

Record 

Variable 

Format 

Default  Purpose 

1 

XQDE 

12 

0  use  default;  other  values 

not  relevant  with  EN967S 

KXOUT 

2X,I2 

0  =  1 ,  print  influence  matrix 

=  0,  otherwise 

Appendix  A 


2c 


Table  A.4 

CONTENTS  OF  EN967D.DAT 

Record  Variable  Format  Default  Purpose 

IH  14  *  number  of  blocks 

NH(I),I-l,f1  1914/2014  *  number  of  control  points 

In  Ith  block 

2  NSB  14  8888  =  0,  full  output  for  last  iteration 

=  8888,  minimum  output 

NOTES 

*  The  defaults  are  NH(I)  »  NSC  and  fl  -  NXC  if  NSC  is  greater  than  20.  Otherwise,  NH(I) 
*  2*NSC  and  n  is  adjusted  accordingly  (NH  (PI)  *  NSC  if  M  is  odd). 


Table  A.5 


CONTENTS  OF  EN967E.DAT 


Record 

Variable 

Format 

Default 

Purpose 

1 

IIBIN 

2X.I2 

0 

»  0,  file  NAME30.DAT  is 
saved 

(see  Table  A.1 ) 
otherwise  not 

1 1  LPT 

12 

0 

=  0,  printer  output  Is  saved 
in  file  EEOUT.DAT, 
otherwise  it  is  printed 
immediately 

SU 

2X.F8.0 

length*breadth 

reference  area  for  force  and 
moment  coefficients,  as  are  B 
and  C 

B 

F8.0 

breadth 

lateral  reference  length 

C 

F8.0 

length 

longitudinal  reference  length 

X.Y.Z 

3F8.0 

0,0,0 

moment  reference  point 

2 

N 

4X,  12 

0 

use  default;  other  values  not 
relevant  with  EN967S 

3 


KDSUB I 


4X.I2 


0 


use  default;  other  values  not 
relevant  with  EN967S 


CPU  TIME  (s) 


Appendix  A 


Figure  A.2.  CPU  Times  for  EN967C/D/E:  DREA  DEC-20  Computer 
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Appendix  B 

USER'S  GUIDE  FOR  VELSPL 


VELSPL  calculates  a  smooth  parametric  B-spllne  approximation  to  the  discrete 
cartesian  flow  vectors  provided  by  EN967.  It  has  been  designed  to  be  used  in 

conjunction  with  the  hull  splining  program  SMHULL1,  the  suite  of  programs  EN967,  and  the 

hull  display  programs  STMHLC  and  DSPHLL10.  VELSPL  Is  equipped  with  a  wide  range  of 
defaults  allowing  both  ease  of  use  for  the  inexperienced  user,  and  generality  for  the 
experienced  user  with  a  difficult  problem. 

VELSPL  is  written  in  FORTRAN  7.0  on  the  DEC-20  at  DREA  and  calls  subroutines  from 
the  spline  libraries  BSPLIN:BSPLN73  and  BSPLIN:HLLYSP4,  from  the  GETWRD  Package11, 

and  from  the  PLOTIO  library12.  At  DREA,  VELSPL  may  be  obtained  by  the  following  calling 
sequence: 

EXECUTE  PS: <HALLY>VELSPL,  PS: <HALLY>SMSUBS/LIB,  PS: <HALLY>GETURD, 
BSPLIN:HLLYSP/LIB,  BSPLIN:BSPLN7/LIB.  PS: <LIB>PL0T18/LIB 

VELSPL  is  designed  to  be  run  on  a  Tektronix  graphics  terminal  or  equivalent, 
although,  if  none  of  the  plotting  features  are  being  used,  it  is  possible  to  run  it  on  a  non¬ 
graphics  terminal. 

VELSPL  interacts  with  the  user  via  a  command  language  described  in  the  next 
section.  The  commands  and  their  functions  are  described  in  subsequent  sections. 


B.1  TERMINALS  AND  PLOTTERS 

In  this  section  information  about  terminals  and  plotters  which  can  be  used  with 
VELSPL  is  presented. 


B.1.1  IDENTIFYING  TERMINAL  AND  PLOTTER  TYPE 

VELSPL  has  been  designed  for  use  with  both  a  terminal  and  a  plotter,  though  in  the 
case  of  PLOTIO  compatible  terminals  these  may  be  one  and  the  same,  and  if  none  of  the 
plotting  features  are  used  the  terminal  need  not  support  graphics.  The  terminal  is  used 
for  all  interaction  between  user  and  program  via  the  command  language,  while  the  plotter 
is  used  for  all  graphic  output.  Immediately  upon  entering  any  of  the  programs,  the  user 
will  be  prompted  for  the  type  of  terminal  being  used.  The  list  below  gives  the  terminals 
currently  recognized  and  the  commands  used  to  identify  them: 

1)  Tektronix  4006:  4006 

2)  Tektronix  4010:  4010 


3)  Tektronix  4012:  4012 
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4)  Tektronix  4013:  4013 

5)  Tektronix  4014:  4014 

6)  Tektronix  401 5:  4015 

7)  Tektronix  401 6:  4016 

8)  Retrographics  upgrade  of  a  Digital  VT100:  VT100/RG 

When  the  terminal  type  has  been  identified,  the  user  is  prompted  to  determine  what 
device  is  being  used  to  generate  the  graphic  displays.  This  may  be  the  terminal  or  a 
plotter.  If  a  plotter  is  used,  the  command  prompts  will  appear  on  the  terminal  screen  but 
all  graphic  displays  are  drawn  on  the  plotter.  Allowed  responses  to  the  prompt  for  the 
type  of  plotter  are 

4662,  4662/31,  4663  and  TERM  I NAL 

where  4662  and  4663  indicate  the  Tektronix  plotters  of  those  numbers,  4662/31  indicates 
a  Tektronix  4662  plotter  with  option  31,  and  TERMINAL  indicates  that  the  terminal  is  being 
used  for  all  graphic  displays.  If  a  plotter  is  being  used,  the  user  will  also  prompt  for  the 
dimensions  of  the  paper  on  which  the  plot  is  to  be  drawn.  These  dimensions  should  to  be 
given  in  inches,  the  horizontal  dimension  first  followed  by  the  vertical  dimension.  An 
optional  command  which  may  be  entered  at  the  samp  time  as  the  page  dimensions  is  the 
extent  of  the  page  to  be  used.  There  are  three  possibilities 

BOTTOM,  FULL  or  TOP 

according  to  whether  the  full  page  is  to  be  used,  or  the  top  or  bottom  half. 

Following  the  prompt  for  the  type  of  plotter,  the  user  will  be  prompted  for  the  baud 
rate  currently  being  used.  Note  that  this  does  not  change  the  baud  rate,  it  only  tells  the 
program  what  the  baud  rate  is. 

A  means  for  changing  the  type  of  terminal  or  plotter  at  any  time  during  the 
execution  of  the  program  is  described  in  Appendix  B.5.7. 


B.1.2  SET  UP  OF  4662  AND  4663  PLOTTERS 

When  a  Tektronix  4662  or  4663  plotter  is  being  used  it  should  be  joined  in  series 
between  the  computer  and  the  terminal.  The  plotter  should  NOT  be  in  Copy  Mode. 
Appropriate  settings  for  the  RS-232-C  switches  on  the  rear  of  the  plotter  are: 

A  »  1 ,  B  =  2,  C  «  2,  D  =  3  for  operation  at  1 200  baud 

A  »  1 ,  B  -  2,  C  =  2,  D  -  1  for  operation  at  300  baud 

Make  sure  the  plotter  is  on  before  starting  the  program.  It  will  not  work  correctly  if  the 
plotter  is  turned  on  Just  prior  to  drawing  the  graphic  display. 


Before  drawing  a  graphic  display,  the  user  will  be  prompted  to  ensure  that  the 
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paper  is  loaded  on  the  plotter.  At  this  time  the  user  should  place  a  clean  sheet  of  paper 
on  the  plotter  having  the  dimensions  specified  by  him  earlier.  The  lower  left  and  upper 
right  plotting  points  should  be  set  at  the  corresponding  points  on  the  paper.  The  program 
will  leave  margins  around  the  display  so  there  is  no  need  to  reduce  the  plotting  area  to 
allow  for  margins.  If  the  plotter  has  Option  31  there  is  no  need  to  load  a  pen:  pen  #1  will 
be  loaded  automatically  by  the  program.  If  it  does  not,  the  user  should  ensure  that  a  pen 
is  properly  loaded  before  proceeding. 


B.2  COMMAND  LANGUAGE  IN  VELSPL 

VELSPL  interacts  with  the  user  via  a  command  language,  implemented  with  the 

GETWRD  package  of  subroutines11.  The  GETWRD  Package  allows  a  type  ahead  facility, 
abbreviated  commands,  command  completion,  help  messages,  and  explicit  error  messages, 
providing  a  friendly  interface  between  the  user  and  the  program.  All  the  special 
character  defaults  in  GETWRD  have  been  chosen  so  that 

a)  The  delimiting  character  is  a  blank.  All  commands  must  be  separated  by 
blanks. 

b)  The  help  character  is  a  question  mark.  For  a  help  message  and  a  list  of 
possible  commands,  type  a  question  mark. 

c)  The  completion  character  is  a  dollar  sign.  To  have  the  current  command 
completed  from  an  abbreviated  form,  type  a  dollar  sign. 

d)  The  delete  letter  character  is  a  back  slash,  \.  To  cause  the  previous  letter 
to  be  disregarded,  type  a  backslash. 

e)  The  delete  word  character  is  an  open  angle  bracket,  <.  To  cause  the 
previous  word  to  be  disregarded,  type  an  open  angle  bracket. 

f)  The  abort  character  is  a  caret,  a.  To  jump  back  to  the  highest  level  of 
command,  type  a  caret. 

At  all  levels  of  command  the  prompt  for  further  input  is  an  arrow,  =>.  To  the  left  of 
the  arrow  appears  a  short  indication  of  the  input  desired.  At  the  uppermost  level  of 
command  the  word  'command'  is  to  the  left  of  the  arrow.  To  the  right  of  the  arrow 
appears  the  text  of  the  command  so  far.  For  example,  the  prompt 

baud  rate«>CHANGE  S-ORDER 

asks  the  user  to  input  a  new  value  for  the  order  of  the  splines  in  the  s-direction.  The 
user  is  at  the  third  level  of  command,  the  previous  two  commands  having  been  CHANGE  and 
S-ORDER.  A  new  user  not  understanding  what  is  required  at  this  point  may  type  a  question 
mark  followed  by  a  carriage  return  to  receive  the  more  explicit  help  message 

Enter  the  order  of  the  splines  in  the  s-direction. 

For  further  details  see  the  GETWRD  Package  Manual11. 


The  command  to  exit  from  VELSPL  is 
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command->EXIT 

In  addition,  VELSPL  has  a  sub-command  mode  in  which  single  character  commands 
are  input  via  the  graphic  cross-hairs  on  the  terminal  screen.  The  appearance  of  the 
cross-hairs  is  the  prompt  for  the  input  of  a  command.  Typing  a  question  mark  will  cause  a 
help  message  to  be  written  listing  the  possible  commands.  The  command  'E'  will  exit  to 
the  higher  command  level.  Further  details  are  given  below. 


B.3  INITIAL  INPUT 

Upon  entering  VELSPL  the  message 
Entering  VELSPL  ... 

is  displayed  followed  by  prompts  for  the  terminal  type,  plotter  type,  and  baud  rate  as 
described  in  Appendix  B.1.  These  are  followed  by  a  prompt  for  the  name  of  the  the 
EN967S  output  file  containing  the  derivatives  of  the  hull  coordinates.  This  file  is  called 
NAME17.DAT  by  EN967S,  but  its  name  might  subsequently  have  been  changed  by  the 
user.  If  no  extension  is  given  the  default  extension  is  .DAT.  The  default  file  name, 
obtained  by  typing  a  carriage  return  only,  is  NAME17.DAT. 

After  reading  the  derivatives  of  the  hull  coordinates,  the  user  is  prompted  for  the 
name  of  the  EN967E  output  file  which  contains  the  cartesian  components  of  the  velocity 
and  the  pressure  coefficients  at  the  panel  control  points.  EN967E  calls  this  file 
NAME30.DAT,  which  is  the  VELSPL  default  name.  Errors  occurring  while  reading  data  from 
either  of  these  two  files  will  be  flagged  and  program  execution  terminated. 

After  reading  the  velocity  components  and  pressure  coefficients,  the  VELSPL 
prompts  for  the  file  in  which  the  spline  coefficients  are  to  be  stored.  The  file  name  can 
be  at  most  six  letters  long  with  a  three  letter  extension:  for  example,  VELSPL.DAT.  If  no 
extension  is  given  the  default  extension  is  .DAT. 

VELSPL  now  enters  command  mode.  The  prompt 
command-> 

is  given,  and  the  program  awaits  a  command  from  the  user. 


B.4  SPLINING  WITH  DEFAULTS 

If  none  of  the  defaults  are  to  be  changed,  splining  of  the  data  is  begun  by  giving 
the  command 

command->START 

VELSPL  begins  by  splining  the  x-components  of  the  contravariant  velocity  along  each 
single  station  followed  by  the  s-components  at  each  station,  then  it  splines  the  spline 
coefficients  of  the  x-components,  and  finally  the  spline  coefficients  of  the  s-components 
(see  Section  5.2).  The  splines  of  the  stations  are  numbered  by  the  corresponding  station 
number.  The  splines  of  the  spline  coefficients  are  numbered  sequentially  according  to 
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the  number  of  the  B-spiine  coefficient  being  splined.  These  numbers  are  referred  to  as 
the  numbers  of  the  spline  rows. 

VELSPL  reports  on  its  progress  by  typing  a  number  on  the  terminal  screen  whenever 
a  spline  is  completed.  For  the  splines  of  the  stations,  the  number  printed  is  the  station 
number.  For  the  splines  of  the  spline  coefficients,  the  number  printed  is  the  number  of 
the  spline  coefficient.  When  it  is  completed,  VELSPL  writes  the  spline  coefficients  in  the 
user-specified  output  file  in  the  following  format: 


Record 

Variable 

Format 

Purpose 

1 

KVS 

13 

Order  of  the  splines  in  the  s-direction 

NVS 

13 

Number  of  B-splines  in  the  s-direction 

2 

TVS  (I) 

5F15.7 

Knots  for  splines  in  the  s-direction 
I-l.KVX+NVX 

3 

KVX 

13 

Order  of  the  splines  in  the  s-direction 

NVX 

13 

Number  of  B-splines  in  the  x-direction 

4 

TVX(I) 

5F15.7 

Knots  for  splines  in  the  x-direction 

I -1 , KVX+NVX 

5 

BCOEFV (IQ, IP,K) 

4F1S.7 

B-spline  coefficients 

IQ-1, NVS;  IP-1, NVX;  K-1,2 

The  following  message  is  then  written: 

Fini shed  spl ining. . . 

If  some  of  the  splines  could  not  acheive  the  accuracy  required  of  them  (i.e.  their  chi- 
squares  are  higher  than  the  input  chi-square  -  see  Appendix  B.4.1)  then  the  following 
message  will  also  be  written: 

Accuracies  of  splines  have  been  written  in  VELADC.DAT 

In  the  data  file  VELACC.DAT,  useful  information  about  the  accuracies  of  the  splines  is 
written.  For  each  of  the  splines  of  insufficient  accuracy,  the  required  accuracy  of  the 
spline  and  the  actual  accuracy  of  the  computed  spline  are  displayed. 

The  program  then  returns  to  command  mode.  Velocity  components  for  the  flow 
around  a  different  hull  may  be  splined  by  giving  the  command 

command«>GET -NEW-OATA 

After  which  the  user  will  be  prompted  for  file  names  containing  new  input  data. 
Otherwise,  exit  from  VELSPL  by  giving  the  command 

cominand->EXI  T 


■ , ' .  - 


%  V 


Appendix  B 


33 


B.5  CHANGING  THE  DEFAULTS 

Although  the  defaults  have  been  chosen  with  care  and  will  work  adequately  most  of 
the  time,  it  is  virtually  impossible  to  have  a  set  of  defaults  which  will  work  for  the  flow 
around  every  hull.  Hence,  VELSPL  allows  the  expert  user  to  change  the  default  values. 

The  defaults  may  be  grouped  into  two  broad  classes:  those  which  may  be  different 
for  each  spline  calculated,  and  those  which  must  be  the  same  for  all  splines  in  a  given 
direction,  the  global  defaults.  An  example  of  the  former  is  the  error  associated  with  any 
particular  data  point.  An  example  of  the  latter  is  the  number  of  knots.  The  global 
defaults  must  be  set  before  any  sptining  begins.  The  remaining  defaults  may  be  changed 
as  the  splining  proceeds.  The  default  values  which  may  be  changed  are  described  in 
detail  in  Appendices  B.4.1  to  B.4.8. 

If  it  is  desired  to  change  defaults  as  the  splining  proceeds  the  command 
command->PLOT -EACH-SPLI NE 

must  be  given.  After  each  spline  has  been  calculated,  it  is  plotted  and  a  sub-command 
mode  entered.  Non-global  defaults  may  be  changed  by  appropriate  sub-commands 
described  below.  The  prompt  in  the  sub-command  mode  is 

sub-command»> 

After  changes  have  been  made  to  the  defaults  at  the  sub-command  level,  the  spline  may 
be  recalculated  by  giving  the  command 

sub-command«>RESPLINE 

Only  the  spline  for  the  current  station  or  spline  row  will  be  recalculated.  If  the  spline  is 
satisfactory,  the  spline  of  the  station  or  spline  row  is  obtained  by  giving  the  sub¬ 
command 

sub-command*>CONT I NUE 

Since  there  are  typically  40  or  50  splines  calculated,  it  is  often  tedious  to  change 
defaults  for  every  spline  when,  perhaps,  only  the  splines  at  one  or  two  stations  will  cause 
difficulty.  For  this  reason,  arguments  may  be  given  to  the  command  PLOT-EACH-SPLINE: 

command-»PLQT-EACH-SPLINE  station/ spline-row  number 

where  station/ spline-row  is  one  of 

STATION  orSPLINE-ROU 

and  number  is  the  number  of  the  station  or  spline  row.  The  splines  will  then  be  calculated 
without  plotting  until  the  indicated  station  or  spline  row  is  reached.  All  splines  will  then 
be  plotted  until  plotting  is  suppressed  by  the  sub-command 

eub-command«>SUPPRESS-PLOTS  station! spline-row  number 

This  command  causes  the  suppression  of  plotting  until  the  indicated  station  or  spline  row 
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is  reached.  If  no  arguments  are  given,  VELSPL  will  complete  its  calculations  without 
plotting  further  splines. 

Once  a  non-global  default  is  changed,  the  new  value  usually  remains  as  the  default 
until  changed  again.  There  are  exceptions  to  this  rule,  however,  so  it  is  usually  best  to 
check  the  default  descriptions  below. 


B.5.1  CHANGING  THE  AMOUNT  OF  SMOOTHING 

The  splines  calculated  by  VELSPL  are  smoothing  splines  and  hence,  the  desired 
amount  of  smoothing  must  be  specified;  that  is,  the  splines  will  not  go  exactly  through  the 
data  points  but  will  pass  close  to  them.  The  distance  the  spline  is  allowed  to  be  from  the 
data  points  must  be  specified.  This  is  done  by  establishing  an  appropriate  chi-square  for 
the  spline.  The  chi-square  is  not  a  global  variable  and  may  be  changed  at  both  the 
command  and  sub-command  levels.  At  both  levels  the  appropriate  command  is 

command«>CHANGE  AMOUNT -OF -SMOOTH I NG  method 

where  method  is  one  of 

FOR-ALL-STATIONS,  FOR-ALL-X-SPLlNES,  or  CALCULATED-FROM-DATA 

The  commands  FOR-ALL-STATIONS  and  FOR-ALL-X-SPLINES  both  require  an  argument: 
a  number  equal  to  the  chi-square  for  each  of  the  splines  of  the  stations  or  the  spline 
rows  respectively. 

The  command  CALCULATED-FROM-DATA  indicates  that  the  chi-square  for  each  spline 
is  to  be  calculated  internally  by  the  function  PRERR  in  the  BSPLIN:HLLYSP  library4.  This  is 
the  default.  Once  changed  the  default  for  the  amount  of  smoothing  remains  at  its  new 
value  until  changed  again.  If  the  chi-square  is  calculated  from  the  data  by  PRERR,  all 
subsequent  splines  will  also  call  PRERR  until  the  chi-square  is  changed  explicitly  by  the 
user. 


B.5.2  CHANGING  THE  DATA  POINT  ERRORS 

Each  data  point  has  an  error  associated  with  it:  the  smaller  ths  error,  the  closer  the 
spline  will  try  to  come  to  its  corresponding  data  point.  Note  that  the  errors  are  relative, 
their  values  being  modulated  by  the  value  of  the  chi-square  of  the  spline  (see  the 
previous  section):  that  is,  by  changing  the  error  of  a  data  point  you  car.  make  the  spline 
come  closer  or  farther  away  from  the  data  point  associated  with  it,  but  the  overall  chi- 
square  of  the  spline  will  not  change.  Changing  all  data  point  errors  by  the  same 
multiplicative  factor  (e.g.  doubling  each  of  them)  will  not  affect  the  spline  at  alt.  The  data 
point  errors  are  non-global  parameters  and  may  be  changed  at  both  command  and  sub¬ 
command  levels.  At  the  command  level  the  errors  can  be  changed  only  by  setting  them  all 
equal  to  a  single  value.  The  appropriate  command  Is 

command->CHANGE  MEAN-ERROR  value 

where  value  is  the  value  which  the  errors  are  to  assume.  Default  for  this  value  is  1.  At 
the  sub-command  level  the  command  for  changing  the  errors  is 
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sub-command«>CHANGE  ERRORS  method. 

where  method  is  one  of 

ALL-EQUAL  or  INTERACTIVE-SPECIFICATION 

The  command  ALL-EQUAL  causes  the  value  of  all  errors  to  be  set  to  the  global  default 
value. 

The  command  INTERACTIVE-SPECIFICATION  allows  the  changing  of  the  errors  of 
individual  data  points  via  the  cross-hair  sub-sub-command  mode  (see  Appendix  B.2). 
Upon  giving  the  command  the  data  points  of  the  current  spline  are  replotted  and  the 
cross-hair  appears  as  a  prompt  for  a  single  letter  command.  The  available  commands  are 

1)  E  or  e  :  Exit  from  the  cross-hair  sub-sub-command  mode  to  the  sub¬ 
command  mode. 

2)  N  or  n  :  A  new  value  for  the  error  of  the  data  point  closest  to  the  cross¬ 
hair  may  be  entered.  A  line  is  drawn  from  the  data  point  identified  to  the 
value  typed  in.  The  line  allows  the  user  to  remember  which  error  value  is 
associated  with  which  data  point. 

3)  O  or  o  :  The  old  value  of  the  error  of  the  data  point  closest  to  the  cross¬ 

hair  is  printed  on  the  screen  somewhat  above  or  below  the  cross-hair 
position.  , 

4)  P  or  p  :  Plots  the  data  points  again.  After  several  errors  have  been 
changed,  the  screen  can  get  cluttered  with  old  and  new  error  values.  This 
command  allows  the  screen  to  be  cleared  of  all  these  numbers. 

5)  ?  :  Causes  a  list  of  possible  commands  to  be  written  above  the  plot. 

The  errors  once  changed  remain  fixed  until  the  program  progresses  from  splining  stations 
to  splining  spline  rows.  At  that  time  the  the  errors  are  all  set  equal  to  their  global  default 
value. 


B.5.3  CHANGING  THE  ORDER  OF  THE  SPLINES 

The  order  of  the  splines  of  the  stations  and  the  order  of  the  splines  of  the  spline 
rows  are  global  parameters  which  may  be  changed  by  the  command 

command»>CHANGE  order-type  order 

where  order-type  is  one  of 

S-ORDER  or  X-ORDER 

and  order  is  the  new  order  of  the  appropriate  spline  type.  X-ORDER  refers  to  the  splines 
of  the  stations  (the  splines  In  the  x-dlrection)  while  S-ORDER  refers  to  the  splines  of  the 
spline  rows  (the  splines  in  the  s-direction). 
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The  order  of  the  splines  should  be  changed  before  changing  the  knots  (see 
Appendix  B.5.4)  as  the  order  of  the  spline  is  used  when  calculating  the  knots.  If  the 
knots  are  changed  before  the  orders  are  changed  a  warning  message  will  be  printed.  The 
default  value  for  the  orders  of  the  splines  in  both  directions  is  4. 


B.5.4  CHANGING  THE  KNOTS 

The  knots  for  the  spline  are  global  parameters,  sc  that  they  may  only  be  changed 
before  the  splining  begins.  The  knots  for  the  splines  of  the  stations  and  the  knots  for  the 
splines  of  the  spline  rows  are  specified  separately. 

The  command  to  change  the  knots  is 

command«>CHANGE  knot-type  method 

where  knot-type  is  one  of 

X-KNOTS  or  S-KNOTS 

and  indicates  that  the  knots  for  stations  or  spline  rows,  respectively,  are  to  be  changed. 
method  is  one  of 

EQUALLY-SPACED,  FROfl-FILE,  or  NOT-A-KNOT 

The  command  EQUALLY-SPACED  requires  an  argument,  the  number  of  knots.  Hence, 
the  command 

command->CHANGE  X-KNOTS  EQUALLY-SPACED  15 

indicates  the  the  splines  of  the  spline  rows  (the  splines  in  the  x-dfrection)  are  to  use  1 5 
equally  spaced  knots. 

The  command  FROM-FILE  indicates  that  the  knots  are  to  be  read  from  a  file.  A 
prompt  is  given  for  the  file  from  which  the  knots  are  to  be  read.  The  knots  are  read  by 
an  unformatted  READ  as  follows: 

READ (UNI T,*)NKT, (TCI), I-l.NKT) 

NKT  is  the  number  of  knots  and  the  array  T  contains  the  values  of  the  knots. 

The  command  NOT-A-KNOT  indicates  that  the  knots  are  to  evaluated  by  the  not-a- 

knot  condition3.  This  condition  ensures  that  the  number  of  B-splines  is  exactly  equal  to 
the  number  of  points  being  splined. 

The  default  for  the  knots  for  the  splines  of  the  stations  is  that  they  be  equally 
spaced  with  the  number  of  knots  equal  to  the  average  number  of  data  points  per  station. 
The  default  for  the  knots  for  the  splines  of  the  spline  rows  is  that  they  be  calculated  by 
the  not-a-knot  condition. 
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B.5.5  CHANGING  THE  STIFFNESS  WEIGHTS 

Each  spline  may  be  altered  by  the  specification  of  input  parameters  called  stiffness 
weights.  There  N-k+1  stiffness  weights  (N  is  the  number  of  B-splines  in  the  spline,  and  k 
is  the  order  of  the  spline)  corresponding  to  the  intervals  between  the  knots  tk,...,tN+1.  By 
increasing  the  stiffness  weight  associated  with  a  given  interval,  the  spline  will  become 
"stiffer"  over  that  interval:  that  is,  the  curve  will  tend  to  be  straighter.  If  the  stiffness 
weight  is  decreased  the  spline  will  become  more  pliant  over  the  associated  interval  and 
will  bend  more  easily.  Spline  wiggles  may  be  avoided  successfully  by  making  the 
stiffness  weights  high  where  the  data  is  very  regular  and  low  where  the  data  is  changing 
rapidly.  Like  the  data  point  errors,  the  values  of  the  stiffness  weights  are  relative: 
changing  them  all  by  the  same  multiplicative  factor  will  not  alter  the  spline  at  all.  More 
complete  information  on  the  stiffness  weights  and  the  manner  in  which  they  are 

implemented  in  the  spline  algorithms  is  given  by  Hally4. 

The  stiffness  weights  are  non-globat  parameters  which  can  be  changed  at  both  the 
command  and  sub-command  levels.  At  the  command  level  give  the  command 

command«>CHANGE  ST I FFNESS-UE I GHTS  method 

where  method  is  one  of 

ALL-EQUAL, INTERACTIVE-SPECIFICATION  or CALCULATED-FRDM-DAT A 

The  command  ALL-EQUAL  sets  all  the  stiffness  weights  to  1 .  Since  only  the  relative 
values  of  the  stiffness  weights  are  important,  there  is  no  need  to  allow  the  setting  of  all 
stiffness  weights  to  a  value  other  than  1 . 

The  command  INTERACTIVE-SPECIFICATION  allows  the  changing  of  individual 
stiffness  weights  via  the  cross-hair  sub-sub-command  mode  (see  Appendix  B.2).  Upon 
giving  the  command  the  data  points  of  the  current  spline  are  replotted  and  the  cross-hair 
appears  as  a  prompt  for  a  single  letter  command.  The  position  of  the  knots  for  the  spline 
are  indicated  by  the  row  of  crosses  at  the  bottom  of  the  plot.  There  is  a  stiffness  weight 
associated  with  each  knot  interval.  The  available  commands  are 

1)  E  or  e  :  Exit  from  the  cross-hair  sub-sub-command  mode  to  the  sub¬ 
command  mode. 

2)  N  or  n  :  A  new  value  for  the  stiffness  weight  of  the  knot  interval  identified 
by  the  vertical  cross-hair  may  be  entered.  A  line  is  drawn  from  the  centre 
of  the  knot  interval  to  the  position  of  the  cross-hair.  The  line  allows  the 
user  to  remember  which  value  is  associated  with  which  knot  interval. 

3)  0  or  o  :  The  old  value  of  the  stiffness  weight  of  the  knot  interval 
identified  by  the  vertical  cross-hair  is  printed  at  the  position  of  the  cross¬ 
hair. 

4)  P  or  p  :  Plots  the  data  points  again.  After  several  stiffness  weights  have 
been  changed,  the  screen  can  get  cluttered  with  old  and  new  stiffness 
weight  values.  This  command  allows  the  screen  to  be  cleared  of  all  these 
numbers. 
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5)  ?  :  Causes  a  list  of  possible  commands  to  be  written  above  the  plot. 

The  command  CALCULATED-FROtt-DATA  indicates  that  the  stiffness  weights  for  each 

spline  are  to  be  calculated  by  the  subroutine  WTIBEG  in  the  BSPLIN-.HLLYSP  library4.  This 
is  the  default. 

New  values  for  the  stiffness  weights  remain  as  defaults  until  changed  again  by  the 
user.  If  the  stiffness  weights  are  calculated  from  the  data  by  WTIBEG,  all  subsequent 
splines  will  also  call  WTIBEG  until  the  stiffness  weights  are  changed  explicitly  by  the 
user. 


B.5.6  CHANGING  THE  SMOOTHING  DERIVATIVE 

The  smoothing  splines  try  to  minimize  a  high  order  derivative  of  the  spline4.  The 
order  of  the  derivative  minimized  is  a  non-global  parameter  which  can  be  changed  at 
either  the  command  or  the  sub-command  level.  The  appropriate  command  is  identical  in 
each  case: 

command->CHANGE  SPIOOTHING-DERIVATI VE  smoothing-derivative 

where  smoothing  derivative  is  the  new  value  of  the  derivative.  In  order  that  the  spline  be 
smooth,  the  derivative  should  have  a  value  greater  than  1 .  Default  value  is  2. 

A  new  value  for  the  smoothing  derivative  remains  as  the  default  until  it  is  changed 
again. 


B.5.7  CHANGING  THE  TERMINAL  AND  PLOTTER  TYPE 

As  described  in  Appendix  B.1,  upon  entering  VELSPL  the  user  is  prompted  for  the 
type  of  terminal  and  plotter  being  used.  If  the  user  wishes  to  change  terminal  or  plotter 
types  at  any  time  (typically  the  user  will  first  check  the  plot  on  the  terminal  screen,  alter 
it  to  his  liking,  then  change  the  plotter  type  from  the  terminal  to  a  plotter  in  order  to 
produce  a  hard  copy)  he  may  do  so  by  giving  the  command 

TERMINAL- 

which  will  trigger  the  same  set  of  prompts  and  responses  as  occurred  at  the  beginning  of 
the  program. 


CALCULATING  STREAMLINES  FROM  CONTRAVARIANT  VELOCITY  COMPONENTS 


Let  XCt)  be  a  streamline  parametrized  by  some  parameter  t:  that  is,  for  all  values  of 
t,  X(t)  lies  on  the  streamline.  Then  dX/dt  is  a  vector  parallel  to  the  flow  vector  V.  A 
vector  perpendicular  to  V  is 

v^vV-vV  (C.1) 


whence  the  condition  that  (x(t),s(t))  is  a  streamline  may  be  written 

".^.0 

dt 

But, 

dX  _  ax  dx  dX  ds  _  1  dx  2  ds 
dt  dx  "  dt  +  ds  '  dt  *  dt  +*  "dt” 


Hence,  equation  (C.2)  may  be  written 


dx 

"dt” 


\,2  _  yl - _  q 


ds 

"dt” 


ds  a 

dX  y2 


(C.2) 


(C.3) 


(C.4) 


(C.5) 


If  the  contravariant  velocity  components  at  each  point  on  the  hull  are  known,  this 
equation  is  easily  integrated  to  yield  the  streamlines. 
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